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We consider a Bose-Einstein condensate subject to a rotat- 
ing harmonic potential, in connection with recent experiments 
leading to the formation of vortices. We use the classical 
hydrodynamic approximation to the non-linear Schrodinger 
equation to determine almost analytically the evolution of 
the condensate. We predict that this evolution can exhibit 
dynamical instabilities, for the stirring procedure previously 
demonstrated at ENS and for a new stirring procedure that we 
put forward. These instabilities take place within the range of 
stirring frequency and amplitude for which vortices are pro- 
duced experimentally. They provide therefore an initiating 
mechanism for vortex nucleation. 
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Quantized vortices in superfluid helium II, in particu- 
lar the issue of vortex nucleation in a rotating container, 
have long been the subject of intense work M. With 
the recent production of gaseous Bose-Einstein conden- 
sates [|| the subject has gained a renewed interest. On 
the experimental side, three groups have succeeded in 
obtaining vortices in atomic condensates, with two dif- 
ferent techniques: a phase imprinting technique at JILA 
H and the equivalent of the helium rotating bucket ex- 
periment at ENS H and at MIT Q . At ENS a rotating 
laser beam superposed onto the magnetic trap holding 
the atoms creates a harmonic rotating potential with ad- 
justable anisotropy e and rotation frequency f2. For a well 
chosen range of variation for f2 one or several vortices are 
nucleated, and then detected as holes in the density pro- 
file of the gas after ballistic expansion Q or by a mea- 
surement of the angular momentum of the condensate 
H . A striking feature of the ENS experimental results is 
that, for a very weak anisotropy e, nucleation of vortices 
takes place in a narrow interval of rotation frequencies 
Pmin,^max] around 0.7ux, where uj± is the mean oscil- 
lation frequency of the atoms in the x—y plane, whatever 
the number of atoms or the oscillation frequency u> z along 
z in the experiment 

While the JILA experiment is well understood theoret- 
ically Q the situation is more involved for the ENS exper- 
iment. Several theoretical articles, inspired by the case of 
superfluid helium, have tried to predict the value of the 
lower vortex nucleation frequency f2 m i n from purely ther- 
modynamic arguments ]lO|-|l(| . The proposed values for 
^min are significantly different from the observed value 
of 0.7w_l, or depend on the trap aspect ratio lo z /lo± or on 
the atom number, in contradiction with the observations 
at ENS. Also thermodynamical reasonings are not able 
to predict the upper vortex nucleation frequency Sl max , 



which is also close to 0.7u)± for low anisotropy e. 

In this paper we consider the time dependent problem 
of a condensate subject to a harmonic stirring poten- 
tial. We use the classical hydrodynamic approximation 
to the time dependent Gross-Pitaevskii equation (GPE) , 
an approximation well justified for the ENS parameters 
|l7| . We are then able to reformulate the partial dif- 
ferential hydrodynamic equations in terms of ordinary 
differential equations, which allows an almost analytical 
solution 0|. Our main result is the discovery of dy- 
namical instabilities in the evolution of the condensate 
for a certain range of the rotation frequency and of the 
trap anisotropy. These instabilities will invalidate the 
classical hydrodynamic approximation after some evolu- 
tion time. We have checked with a numerical solution 
of the Gross-Pitaevskii equation that vortices then enter 
the condensate. 

The existence of such a dynamical instability ex- 
plains why in earlier numerical work the time dependent 
Gross-Pitaevskii equation was found to nucleate vortices 
|l4],[l9 20 1 . Furthermore the instability range that we pre- 
dict is very close to the experimentally observed range of 
vortex nucleation, for various stirring procedures. For 
the stirring procedure of ^,^| we recover the "universal" 
numerical value 0.7 for leading to vortex nucleation 

for low anisotropies e. We provide a simple physical inter- 
pretation of this value: for f2/wj_ = l/v^ — 0.7 the har- 
monic stirring potential resonantly excites a quadrupole 
mode of the condensate 0, which induces large oscilla- 
tions of the condensate and eventually a dynamical in- 
stability sets in. We also investigate a new excitation 
procedure to nucleate vortices, that has recently been 
implemented at ENS pi] : the rotation of the stirring 
potential is set up very slowly. The gas then follows 
adiabatically a branch of steady state until the branch 
becomes dynamically unstable. The corresponding lower 
rotation frequency f2 that we predict for vortex nucle- 
ation is also very close to the experimental value. 

In our model atoms are trapped in a harmonic poten- 
tial rotating at the instantaneous frequency Q(t) around 
z axis. For convenience all the calculations of this pa- 
per are done in a rotating frame where the trap axes are 
fixed. The trapping potential then reads: 



u (r,t) = -mcJ^ 



[1 - e(t)]x 2 + [1 + e(t)]y 2 + -L 



(1) 



where m is the mass of an atom, e(t) is the trap 
anisotropy at time t. The parameters lo± and uj z are the 
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oscillation frequencies of the atoms in transverse and ax- 
ial directions for vanishing anisotropy of the stirring po- 
tential. Within the mean field approximation, the time 
evolution of the condensate field or macroscopic wave- 
function ip(r,t) can be described by the time dependent 
Gross-Pitaevskii equation |22[ : 



dip 



-^V 2 + U(r,t) + g\^\ 2 - ft(i)4 
Zm 



^ (2) 



where g = 4irh a/m is the coupling constant, propor- 
tional to the s-wave scattering length a of the atoms, here 
taken to be positive, and where the inertial term propor- 
tional to the angular momentum operator L z along z-axis 
accounts for the frame rotation. The condensate field tj) 
can be written in terms of density p and phase S 7 



(3) 



The equation obtained from the GPE for p is just the 
continuity equation. The equation for S contains the so- 
called quantum pressure term h 2 V 2 y/p/2m^/p that we 
neglect here as compared to the mean-field term pg in 
the Thomas-Fermi approximation. We obtain: 



dp 
di 

dS_ 




n(t) x r 



(4) 



2m 



+ U(r, t)+gp- (ft (i) x f ) • VS. (5) 



A very fortunate feature of the harmonic trap is that 
these superfluid hydrodynamic equations can be solved 
exactly for a condensate initially at equilibrium in the 
non-rotating trap with the followingquadratic ansatz for 
the condensate density and phase |23]: 



Pc(f,t) = p (t) 
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S c (r,t) = s (t) Imui ^ XiBij(t)xj, 



(G) 
(7) 



where Xi, x-i and X3 are the coordinates along x, y and 
z axes respectively. The time dependent dimensionless 
coefficients A^ and By form 3x3 symmetric matrices A 
and B which from Eqs.(Q,^) obey the evolution equations: 



_ 1 dA 



ft 



to , x — = -2A TrB - 2{A, B} H [R, A], 

at u> 1 



^dB 
- ~dt 



-2B 2 -W -A 



O 



\R,B] 



(8) 
(9) 



where {, } stands for the anti-commutator, [, ] stands for 
the commutator of two matrices, the matrix W is di- 
agonal, with components W\\ = (1 — e(i))/2, W22 — 
(1 + e(*))/2, and W33 = (w z /w ± ) 2 /2, and the matrix R, 



originating from the vectorial product in L z , has vanish- 
ing elements except for R12 = —R21 = 1 |24|. Note that 
these equations do not depend on the number of atoms 
nor on the coupling constant g. 

In a first stage, it is very important to study steady 
state solutions of the above equations. We restrict to so- 
lutions that have the same symmetry as the initial state: 
even parity along z, this parity being preserved by time 
evolution. We then find an unique class of solutions, re- 
producing the results of |^5| : the condensate phase varies 
as S(f) — mujj_Pxy where (3 is a real root of 



n 2 



0- 



ft 

UJ ± 



0. 



(10) 



The steady state matrix A is diagonal with elements 
given in [ p5[ . We have plotted in figure |l| the values 
of as function of the rotation frequency ft for a fixed 
anisotropy e. For ft between zero and a bifurcation value 
ftbif(e) depending on e there is a single branch of solu- 
tion for [3. This branch is supplemented by two extra 
branches when ft > ftbif(e). 

We now turn back to the time dependent problem. 
Clearly a condensate with a vortex cannot be described 
within the quadratic ansatz (||J^) as the phase S c corre- 
sponds to an irrotational velocity flow. The actual sce- 
nario for the vortex nucleation that we put forward is the 
following: initially very small deviations 5p(r,t) of the 
condensate density and SS(r, i) of the condensate phase 
from the quadratic shapes p c and S c may grow expo- 
nentially fast in the course of time evolution, eventually 
leading the condensate to a structure very different from 
Eqs.(||,^|). This may happen when a dynamical instabil- 
ity is present. 

To reveal such an instability we obtain from the evo- 
lution equations (|], ^|) linearized equations of motion for 
initially small deviations Sp and SS from p c and S c : 



DSp ,. / VSS\ t V 2 S C 

— — = -div p c - dp 

Vt \ m I m 



DSS 
Dt 



-gSp. 



(11) 
(12) 



In these equations, we have introduced the convective 
jL+v c {?,t)-V where v c = VS c /m— ftxr 



derivative -Si 



is the velocity field of the condensate in the rotating 
frame. A polynomial ansatz for SS and Sp of an arbi- 
trary total degree n in the coordinates x, y and z solves 
these linear equations exactly |26[ . This is another nice 
consequence of the harmonicity of the trap. In practice, 
we calculate the evolution operator U n (t) mapping the 
coefficients of the polynomials at time zero onto their 
values after a time evolution t. Dynamical instability 
takes place when one or several eigenvalues of U n grow 
exponentially fast with time t. Note that after rescaling 
of the variables, Eqs.(fTlJTJ) become independent of the 
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number of atoms and of the coupling constant g, in a way 
similar to Eqs.(|]]|). 
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FIG. 1. Phase parameter /3 for a steady state condensate 
as function of the rotation frequency SI in units of wj_ ■ Dotted 
lines: e = 0. Solid lines: e = 0.022. For e = 0.022 and 
uj z /uj± — 0.1 the thick lines on the curves indicate where the 
solution has a dynamical instability (a) of degree n — 2 and 
(b) of degree n = 3. The vertical line is the border of the 
center of mass instability domain for e = 0.022. 

Now we perform a linear stability analysis for two dif- 
ferent stirring procedures of the condensate. 

Procedure I: The ellipticity of the stirring potential 
e is kept fixed and the rotation frequency Sl(t) of the 
stirrer is very slowly ramped up from zero to its final 
value. The condensate, initially in steady state with a 
vanishing parameter (3, adiabatically follows the upper 
branch of steady states with j3 > 0, see figure It is 
then sufficient to determine the dynamic stability of the 
upper branch of steady states. This greatly simplifies 
the calculation as one just has to identify eigenmodes of 
Eqs.( pd| , ^2[ ) evolving in time as exp(Ai) where the eigen- 
value A is a complex number. Dynamical instability takes 
place when A can have a strictly positive real part. As 
shown in figure [j] for e = 0.022 we find that the upper 
branch for j3 is stable for modes of degree n = 2 but 
presents two instability intervals for the modes of degree 
n = 3, a very narrow interval around SI = 0.58a;j_ and a 
broader interval starting aX SI = 0.778wj_. 

We have investigated in a systematic way the instabil- 
ity range of the upper branch of steady states, by vary- 
ing the anisotropy e, the rotation frequency SI and the 
degree n. The instability domain in the SI — e plane for 
modes of degree 3 is mainly made of a crescent, and the 
inclusion of higher degree modes (n = 4, 5) add extra 
crescents from above, see figure ||a. Each crescent has 
on the e = axis (i) a broad basis at SI > uj±/^/2, with 
a non-zero instability exponent, and (ii) a very narrow 
edge at SI = uj± / s/n < uj± / \[2 with a vanishing insta- 
bility exponent. We show in figure ||b that the maximal 
instability exponent for n — 3 has a remarkably weak 
dependence on uj z /lu±. 




FIG. 2. For the upper branch of steady state condensates: 
(a) Dark areas: instability domain in the Q — e plane for 
lj z /lu± =0.1 for degrees n equal to 3, 4 and 5 (crescents from 
bottom to top). There is no dynamical instability for n < 2. 
Solid line: border SI 2 — (1 — °f the branch existence 

domain, (b) Maximal instability exponent Re(A) for n = 3 as 
function of SI, for e = 0.022, and uj z /lj ± = 0.1 (o), 0.5 (o), 1.0 
(x) and y/8 (+). SI, Re(A) are in units of lo±. 

Procedure II: This is the original experimental sce- 
nario of |4||| , where the stirring potential is rotated at a 
fixed frequency and the ellipticity of the stirrer is turned 
on from zero to its final value e/ abruptly. In this case 
we cannot rely on adiabatic following for the condensate 
density and phase, so that we solve the time dependent 
equations for p c (t) and S c (t). Then we perform 

a linear stability analysis as discussed above: we evolve 
a generic polynomial ansatz of degree n for the fluctua- 
tions Sp and 8S according to Eqs.( |n|Jl2] ), which allows to 
construct the evolution operator M n (t) and to calculate 
Zmaxit), the eigenvalue of U n {t) with the largest modu- 
lus. Then we define the mean instability exponent Re (A) 
as the mean slope of In |Z max (£)l as function of time. 

This reveals that within certain range of rotation fre- 
quency the system becomes dynamically unstable, see the 
solid line in figure |^. In the limit of a low anisotropy e the 
instability sets in when the rotation frequency SI is close 
to the value ~ 0.7uj±: in the lab frame, the stirring poten- 
tial of frequency 2SI is then resonant with a quadrupole 
mode ]i~7| , p7| of the condensate of frequency \/2^±, and 
induces large amplitude oscillations of the condensate, 
resulting in a dynamical instability. More precisely, the 
condensate described by the quadratic ansatz p c , S c has 
an angular momentum oscillating in time around a non- 
zero value L z . This value L z is a peaked function of the 
rotation frequency SI, shown in dashed line in figure 3. 
The peak of L z is not exactly located at SI = 0.7oj± be- 
cause of non- linear effects in Eqs.(||, ||). The peak struc- 
ture of the instability exponent in figure 3 is alike the 
peak structure of L z , with a narrower width as dynami- 
cal instability of the vortex free solution p c , S c sets in for 
the higher values of L z only. For values of SI significantly 
above or below 0.7wj_ the stirrer is out of resonance with 
the quadrupole mode and induces only small and stable 
oscillations of the condensate. For larger values of e, the 
instability interval in SI broadens. We have also checked 
that the instability interval depends weakly on uj z /lo±. 
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FIG. 3. For the stirring procedure II, with e/ = 0.01 and 
lj z /uj± = 0.1: Solid line: mean instability exponent Re (A) 
(see text) of the vortex free classical hydrodynamic solution 
p c , S c as function of £7, for n = 3. Dashed line: mean angu- 
lar momentum per particle L z (see text) obtained from p c , S c . 
Filled disks: experimentally measured angular momentum L z 
per particle in the condensate after vortices have possibly en- 
tered the condensate J8j. The initial steady state condensate 
in the calculation of L z has a chemical potential [i = 10hui±, 
close to the experimental value. Re (A) and Q, are in units of 
u±, and L z , L z are in units of Ti. 

What is the connection between the dynamical insta- 
bilities found here and the nucleation of vortices ? To 
obtain a theoretical answer to this question, one has to 
go beyond a linear stability analysis to determine the 
evolution of the condensate in the long run: for a few 
values of the rotation frequency Q and for procedures 
I and II, we have checked by a numerical integration 
of the time dependent GPE in three dimensions, that 
vortices are indeed nucleated in the predicted instability 
domains: after some evolution time, the angular momen- 
tum in the numerical solution suddenly becomes larger 
than the classical hydrodynamic prediction, as vortices 
enter in the condensate. An experimental answer to this 
question for the stirring procedure I has been provided 
recently at ENS | pl[ : the first clear evidence of a vor- 
tex appears for a rotation frequency £1 = 0.77u>_l, very 
close to our prediction ft = 0.778^^. The agreement is 
also very good for procedure II as shown in figure [5| our 
instability domain in Q coincides with the experimental 
vortex nucleation interval within a few percent. 

In summary, the dynamical instabilities that we have 
identified provide an initiating mechanism for the pro- 
duction of vortices in a condensate stirred by a harmonic 
potential, in excellent agreement with the experimental 
results at ENS. 
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